#核函数
import time
t=time.time()
def kernel(a,b,kernel_type=None):
    if kernel_type==None:
        return numpy.dot(a,b.T)
def clipAlpha(alpha,H,L):
    if alpha >=H:
        alpha = H
    if L >= alpha:
        alpha = L
    return alpha
def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(numpy.random.uniform(0,m))
    return j
from matplotlib import pyplot
import numpy
numpy.random.seed(100)
data=numpy.loadtxt("./testSet.txt",delimiter="\t")
# pyplot.scatter(data[:,0],data[:,1],c=data[:,-1],edgecolors="red",linewidths=1.5,alpha=0.7)
# pyplot.show()
label=data[:,-1]
feature=data[:,:-1]
m,n=feature.shape
#给定罚因子
C=100
#初始化alpha和b
b=0
alpha=numpy.zeros(m)
#查找精度
tol=1e-3
#最大迭代次数
max_iters=40
#计算初始误差
E=numpy.sum(label*alpha*kernel(feature,feature),axis=1)+b-label
#开始选择优化变量,检查是否满足KKT条件
iters=0
c=0
while iters<max_iters:
    changed=0
    for i in range(m):#扫描整个数据集
        #第一个变量的选择
        if ((label[i] * E[i] < -tol) and (alpha[i] < C)) or ((label[i] * E[i] > tol) and (alpha[i] > 0)):
            #第二个变量的选择
            #j=numpy.argmax(numpy.abs(E-E[i]))#仅仅是这样不收敛
            # if j==i:
            #     continue
            # elif E[j]==0:
            #     j = selectJrand(i, m)
            #print(E[j]-E[i])
            #计算新的alpha的值
            j = selectJrand(i, m)
            alphaj_old = alpha[j].copy()
            alphai_old = alpha[i].copy()
            eta=kernel(feature[i],feature[i])+kernel(feature[j],feature[j])-2*kernel(feature[i],feature[j])
            #print((E[i] - E[j]) / eta)

            if eta<=0:
                continue
            alphaj_new=alpha[j]+label[j]*(E[i]-E[j])/eta

            if (label[i] != label[j]):
                L = max(0, alphaj_old - alphai_old)
                H = min(C, C + alphaj_old - alphai_old)
            else:
                L = max(0, alphaj_old + alphai_old - C)
                H = min(C, alphaj_old + alphai_old)
            if L==H:
                continue
            alpha[j] = clipAlpha(alphaj_new, H, L)
            #alpha[j]=alphaj_new

            alpha[i]=alphai_old+label[j]*label[i] * (alphaj_old - alpha[j])
            b1=b-E[i]-label[i]*kernel(feature[i],feature[i])*(alpha[i]-alphai_old)-label[j]*kernel(feature[i],feature[j])*(alpha[j]-alphaj_old)
            b2=b-E[j]-label[i]*kernel(feature[i],feature[j])*(alpha[i]-alphai_old)-label[j]*kernel(feature[j],feature[j])*(alpha[j]-alphaj_old)
            if (0 < alpha[i]) and (C > alpha[i]):
                b=b1
            elif (0 < alpha[j]) and (C > alpha[j]):
                b=b2
            else:
                b=(b1+b2)/2
            E = numpy.sum(label * alpha * kernel(feature, feature), axis=1) + b - label
            changed=changed+1
            print(b)
            c=c+1
    if changed==0:
        iters=iters+1
print(b)
print(time.time()-t)
print("OK")
pyplot.scatter(feature[:,0],feature[:,1],c=label)
pyplot.scatter(feature[:,0][alpha>0],feature[:,1][alpha>0],c="none",edgecolors="red",linewidths=1.5,alpha=0.7,s=150)
pyplot.show()
#优化的SMO算法